Correlated electron states and transport in triangular arrays 
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We study correlated electron states in frustrated geometry of a triangular lattice. The interplay 
of long range interactions and finite residual entropy of a classical system gives rise to unusual 
effects in equilibrium ordering as well as in transport. A novel correlated fiuid phase is identified 
in a wide range of densities and temperatures above freezing into commensurate solid phases. The 
charge dynamics in the correlated phase is described in terms of a height field, its fluctuations, and 
topological defects. We demonstrate that the height field fluctuations give rise to a "free" charge 
flow and finite dc conductivity. We show that freezing into the solid phase, controlled by the long 
range interactions, manifests itself in singularities of transport properties. 
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I. INTRODUCTION 

The properties of geometrically frustrated systems are 
complex due to the presence of an energy landscape with 
many degenerate or nearly degenerate minima.^ These 
systems exhibit qualitatively new effects, such as an ex- 
tensive ground state degeneracy and suppression of freez- 
ing down to zero temperature. The latter phenomena 
traditionally have been studied in the context of antifer- 
romagnetic spin models on appropriately chosen lattices, 
and other models with short range interactions. 

Below we demonstrate that similar effects of geomet- 
rical frustration can naturally arise in a classical charge 
system. Namely, we consider ordering and dynamics of 
classical electrons on a triangular lattice with repulsive 
Coulomb interaction between charges on different lattice 
sites. The frustration reveals itself in the phase diagram, 
and also makes charge dynamics and transport very un- 
usual. Here we focus on the electric transport as a natural 
means to study ordering types and phase transitions in 
a charge system. 

Two-dimensional frustrated lattices arise in a variety 
of experimental systems. Artificial structures, such as 
Josephson junction arrays^ and arrays of quantum dots,^ 
have recently become available. One attractive feature of 
these systems is the control of the Hamiltonian and, in 
particular, of frustration, by the system design. Also, ex- 
perimental techniques available for probing magnetic flux 
or charge ordering, such as electric transport measure- 
ments and scanning probes, are more diverse and flexi- 
ble than those conventionally used to study magnetic or 
structural ordering in solids. There have been extensive 
theoretical'''^ and experimental^ studies of phase transi- 
tions and collective phenomena in Josephson arrays. An- 
other example of frustrated charge systems is provided by 
novel superconducting materials based on Co02.^ 

In the present work we discuss a realization of a frus- 
trated charge system in a 2d triangular array of quan- 
tum dots. Recent progress in the epitaxial and litho- 
graphic techniques made it possible to produce^ regu- 



lar and irregular arrays of quantum dots, in which the 
size of an individual dot can be tuned in the 10-100 nm 
range with the rms size distribution of 10-20%. Such ar- 
rays are made of InAs or Gc islands embedded into the 
semiconductor substrate. This stimulated experimental^ 
and theoretical* investigation of electronic ordering and 
transport in these arrays. Capacitance and conductivity 
measurements^ show the quantized nature of charging 
of the dots. However, the interdot Coulomb interaction 
in such systems is weak compared to the individual dot 
charging energy as well as potential fluctuations due to 
disorder in the substrate. 

A promising system fabricated and studied recently^ 
involves nanocrystallite quantum dots which are synthe- 
sized with high reproducibility, of diameters ~ 1.5 — 15 
nm tunable during synthesis, with a narrow size distribu- 
tion (< 5% rms). These dots can be forced to assemble 
into ordered 3d closely packed colloidal crystals,^ with 
the structure of stacked 2d triangular lattices. High flex- 
ibility and structural control open a possibility to study 
effects inaccessible in the more traditional self-assembled 
quantum dot arrays fabricated using epitaxial growth 
techniques. In particular, the high charging energy that 
can reach or exceed the room temperature scale, and the 
triangular lattice geometry of the dot arrays^ are of in- 
terest from the point of view of exploring novel aspects 
of charge ordering and transport. 

From a theoretical viewpoint, charge ordering is closely 
related to, or can be interpreted in terms of, a suit- 
ably chosen spin system. Spin models, which serve as a 
paradigm in a theory of critical phenomena, have appli- 
cations to ordering in different systems, such as the phase 
transitions in adsorbed monolayers. In our analysis we 
map the charge problem onto the triangular antiferro- 
magnetic Ising spin problem. In the situation of interest, 
when charging energy enforces single or zero occupancy 
of the sites, one can interpret the occupied and unoccu- 
pied dots as an 'up spin' and 'down spin' states. Since 
the like charges repel, the corresponding effective spin 
interaction is indeed of an antiferromagnetic kind. Ap- 
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propriately, charge density plays the role of spin density, 
and the gate voltage corresponds to an external magnetic 
field. Besides, one can map the offset charge disorder 
(random potentials on the dots) onto a random magnetic 
field in the spin problem. 

There are, however, a few theoretical and experimen- 
tal aspects of the charge-spin mapping that make the two 
problems not entirely equivalent. First, charge conserva- 
tion in the electron problem gives rise to a constraint on 
total spin in the associated spin problem. This leads to a 
dynamical constraint, namely blocking of single spin flips. 
Microscopically, spin conservation requires the Kawasaki 
(or type B) dynamics^'^ as opposed to the nonconserving 
Glauber (or type A) dynamics. This makes no differ- 
ence with regard to the thermodynamic state at equilib- 
rium, since the system with fixed total spin is statistically 
equivalent to the grand canonical ensemble. However, the 
order parameter conservation manifests itself both in a 
slower dynamics^^ and in collective transport properties 
of the correlated fluid phase discussed below. 

Another important difference between the charge and 
spin problems is in the form of interaction. Spin systems 
are visually desciribed by a nearest neighbor interaction. 
In the spin problem relevant for this work, the triangular 
Ising antiferromagnet (AIAFM ) with the nearest neigh- 
bor interaction, an exact solution in zero field has been 
obtained by Wannier^^ who demonstrated that there is 
no ordering phase transition at any finite temperature. 
In contrast, in the charge problem studied in this arti- 
cle the long range Coulomb interaction makes the phase 
diagram more rich. We find phase transitions at finite 
temperature for certain charge densities. 

One of the main objects in our focus is the corre- 
lated fluid phase, that arises at relatively low tempera- 
tures due to geometric frustration preventing freezing. In 
this phase, equilibrium fluctuations and transport exhibit 
strong correlations. We employ a nonlocal description in 
terms of the height field, used earlier in spin models, 
to map charge dynamics onto Gaussian fluctuations of 
the height surface in the presence of topological defects 
(dislocations). We find that both the height field fluc- 
tuations and the presence of defects contribute to low- 
temperature transport properties. 

The outline of the paper is as follows. In Section II a 
charge Hamiltonian is introduced and stochastic Monte 
Carlo (MC) dynamics is defined, based on the correspon- 
dence between the charges and Ising spins. This dynam- 
ics is used to obtain the phase diagram (Section III) and 
to study the dc conductivity as a function of temperature 
and electron density (Section IV). 

In Section V a height field order parameter is de- 
fined and used to describe the charge ordering. In Sec- 
tion VI we compare the contributions to transport due 
to the height surface fluctuations and the topological de- 
fects. Subsequently, in Section VII, we study the height 
variable fluctuations and evaluate their effective stiffness 
from the MC dynamics. By comparing the stiffness to 
the universal value we rule out the Berezinskii-Kosterlitz- 



Thouless transition and prove that the defects are al- 
ways unbound. In Section VIII we present scaling argu- 
ments that support the conclusions of Section VII, and 
set bounds on the height surface stiffness. In Section IX 
we consider dynamical fluctuations of the height variable 
and estimate their effect on dc conductivity. 

II. THE MODEL 

The Hamiltonian TieX of a quantum dot array describes 

charges qi on the dots, their Coulomb interaction as well 
as coupling to the background disorder potential (/>(r) and 
to the gate potential V^: 

i,j Ti 

Here the vectors run over a triangular lattice with the 
lattice constant a, and Vij = Vi — Vj. The interaction 
V{r) includes a term describing screening by the gate: 

Here e is the dielectric constant of the substrate, and d/2 
is the distance to the gate plane. (The parameter d in 
(2) which controls the interaction range is chosen in our 
simulation in the interval < d/2 < 5a.) The single dot 
charging energy 5V^(0) = e^/2C is assumed to be large 
enough to inhibit multiple occupancy, so that = 0, 1. 

The exponential factor in (2) is introduced for conve- 
nience, to control convergence of the sum in (1). Below 
we use = 2d. In the case of spatially varying e the 
form of interaction can be more complicated. For exam- 
ple, for an array of dots on a dielectric substrate, one 
should replace e in Eq. (2) by (e -I- l)/2. 

Electron tunneling in the system^"^^ presumably oc- 
curs mainly between neighboring dots. The tunneling is 
probably incoherent, i.e. assisted by some energy relax- 
ation mechanism, such as phonons. Since the tunneling 
coupling of the dots is weak,^^^ we focus on the charge 
states and ignore the effects of electron spin, such as ex- 
change, spin ordering, etc. 

The incoherent nature of electron hopping warrants 
employing stochastic MC dynamics in which transitions 
between different states take place in accordance with 
Boltzmann probabilities. The states undergoing the MC 
dynamics are charge conflgurations with = 0, 1 on a 
N X N patch of a triangular array. Periodic boundary 
conditions are imposed by extending the N x N config- 
urations Qi along with the Hamiltonian (1) periodically 
over the entire plane. Spatial periodicity of the MC dy- 
namics is maintained by allowing charge hops across the 
boundary, so that the charges disappearing on one side 
of the patch reappear on the opposite side. 

Charge conservation gives a constraint '^Qi = const 
that has to be enforced throughout the MC evolution. 
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While contributing to MC dynamics slowing down, this 
constraint has no effect on the statistical equilibrium and 
the thermodynamic properties. To study ordering one 
can employ the nonconserving A dynamics, with the the 
gate voltage Vg serving as a parameter controlling the 
system state. This is beneficial due to relatively high 
speed of the A dynamics. In turn, the somewhat slower 
charge conserving B dynamics is used to investigate con- 
ductivity at fixed charge density 

n = A-^J2qi, A = N\ (3) 

i 

Below we recall the definition of the dynamics A and 
B, and introduce the charge-spin mapping which will be 
used throughout the paper. 

A. The dynamics of type A 

In the MC dynamics of type A, since charge is not 
conserved, the charge state is updated independently on 
different sites. The occupancy of a randomly selected site 
i is changed or preserved with the probabilities Wj and 
Wi, which depend on the system state and its attempted 
change as follows: 

Wi/Wi = eM-5E\^^ iT.x) . (4) 
8e\^'^ = (5) 

where Wi + Wi = 1. Here Sqi is the attempted change of 
the charge at the site i, whereby the potential is 

^i= J2v{vij)qj+V^ + cP{vi), (6) 

and Tel is the temperature in the charge system (for 
brevity, from now on we set fcs = !)■ At each MC time 
step, the charge configuration is affected only on one site. 
The potential $i at this site is obtained using the current 
system state modified at the preceding MC step. 

B. The dynamics of type B 

To enforce charge conservation appropriate for the B 
dynamics, we update the state of the system by exchang- 
ing charges on randomly chosen neighboring sites. Specif- 
ically, we randomly select a site i. Then, at each MC time 
step, we randomly chose a site j neighboring to i until the 
occupancy of the sites i and j are not the same, qj ^ qi. 
After that, the occupancies of the sites i and j are ex- 
changed with probability Wi^j, and remain unchanged 
with probability Wi^i, such that 

Wi^j /Wi^i = expi-SE^^p /Tel) , (7) 
6eIP = (g.-g,)($,-$.)-(g.-5,)2y(r,.), (8) 

where Wi^j + Wi^i = 1, and the potential $i is defined 
by Eq. (6). 



C. Charge-spin mapping 

The model (1),(2),(4),(7) possesses an electron-hole 
symmetry, exhibited by introducing a 'spin' variable 

Si = 2qi-l = ±l. (9) 

The charge Hamiltonian (1), rewritten in terms of the 
spin variables (9), becomes 

Wei = + const, (10) 

^« = lT.vi^i3>i^j+T.if^+^^i''i)>i- (11) 

Here we introduced the chemical potential 

^ = 214 + Vu^o (12) 

that corresponds to an external field for the spins s,. 
Here Vk = J2j e^'^^^V{rj) is the Fourier transform of the 
interaction (2). In terms of spin variables the charge 
density (3) is given by 

n = A-'Y.^{si + l), A = N\ (13) 

i 

In the corresponding stochastic dynamics for the spin 
system, defined as above, we use Si instead of qi in the 
energies (5) and (8), and rescale Tei to the spin temper- 
ature 

T = 4Tei. (14) 

Eqs. (9) - (14) provide a mapping of the charge prob- 
lem onto a spin system with long-range interaction (2). 
The positive sign of the coupling (2) corresponds to anti- 
ferromagnetic nearest neighbor interaction. We note that 
the limit d <C a corresponds to the AIAFM problem^^"^^ 
with nearest neighbor interaction. Unless explicitly 
stated, below we consider a system without disorder, 
</.(r) = 0. 

III. PHASE DIAGRAM 

A. Cooling curves 

To explore the ground states as a function of the chem- 
ical potential /x, we use the type A (nonconserving) MC 
dynamics. To reach the true equilibrium at low temper- 
atures the usual precautions are taken by running MC 
first at an elevated temperature, and then gradually de- 
creasing it to the desired value. 

The cooling curves in Fig. 1 show the temperature de- 
pendence of electron density n{T) for a moderate value 
of the screening length d = 2a, with the temperature T 
measured in imits of the nearest neighbor coupling V{a). 
Due to the electron-hole symmetry n <-^- 1 — n, it suffices 
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Temperature 



FIG. 1: (Color online) Cooling curves n{T) at fixed gate volt- 
age for the screening length parameter d — 2a. Due to the 
electron- hole symmetry, only the densities < n < 1/2 are 
shown. The values of fi, related to Vg via (12), are given on 
the right side of the plot. Not labelled are the curves with 
fi = 0.1, 0.2, 0.3 converging to n = 1/2. Additional cooldowns 
are shown for /i = 2.5 and fi = 3.0 (dashed line). 



to consider only < n < 1/2. Similar curves obtained 
for the AIAFM problem realized at small d ^ a are dis- 
played in Fig. 2. 

The two families of cooling curves are qualitatively 
similar in the character of temperature dependence: slow 
at high T, faster at lower T, and exhibiting strong fluc- 
tuations before final stabilization at the T = value. 
Also, in both cases the trajectories are attracted to the 
densities n — and n = 1/3, 2/3. We note, how- 
ever, qualitatively new features in the case of long range 
interaction (Fig. 1). The most obvious one is that the 
values of n attained at T — > span a range of intermedi- 
ate densities, with n being a continuous function of Vj,. 
This behavior, indicating freezing transitions at all den- 
sity values, should be contrasted with the T — > behav- 
ior in the AIAFM case, characterized by discontinuous 
jumps between n = 0, 1/3, 2/3, 1. The latter four values 
correspond to the incompressible states [plateaus in the 
dependence of n vs. Vg] . 

The long range interactions give rise to ordering at 
the densities not realized in the AIAFM problem, the 
simple fraction n = 1/2 being the most prominent one 
(Fig. 1). This density is an attractor for a family of 
cooling curves at ^ 0.5: Several such curves, obtained 
for /i — 0.05, 0.1, 0.2, 0.3, are shown in the top part of 
Fig. 1. Evidently, the basin of attraction of the n = 1/2 
state is considerably smaller than that of the n = 1/3 
and n = 2/3 states, which is to be expected, since the 
ordering at n = 1/2 is controlled by the next-to-nearest 
neighbor interactions. 




Temperature 



FIG. 2; (Color online) A family of cooling curves obtained in 
the same way as in Fig. 1 for d <^ a (the pure AIAFM case) . 
The values of /i are given on the right. 

The result of coohng, in general, is found to depend 
somewhat on the cooling history, especially near the in- 
compressible densities. For different initial random dis- 
tributions of charge, and depending on the specific se- 
quence of MC moves, MC relaxation can lead to differ- 
ent ground states. This happens because for a generic 
long-range interaction there are many states nearly de- 
generate in energy, as illustrated in Fig. 1 by the two 
pairs of curves starting at /i = 2.5 and /i = 3.0 obtained 
for different runs of the MC dynamics. These curves, 
identical at high T, diverge below the temperature inter- 
val where fluctuations develop. The dependence on the 
cooling history in the fluctuation region limits the accu- 
racy of phase diagram obtained from cooling trajectories 
in a finite system. 

Figure 3 summarizes in a schematic way the results of 
the MC study of cooling. It shows the phase diagram 
of the system in the (/i, T) plane. Due to the particle- 
hole equivalence, it posesses the /i ^ — /i symmetry. For 
a generic value of the chemical potential we find three 
distinct temperature phases: the disordered state at high 
temperature, the correlated fluid phase at intermediate 
temperatures, and the solid phases at low T. 



B. Freezing transitions 

Ordering at different n resembles, and indeed can 
be connected to, the phase transitions in adsorbed 
monolayers. The latter have been mapped on the 
known statistical models (Ising, Potts, etc), some of 
which are exactly solvable. Possible phase transitions 
in 2d have been classified by Domany et aL^° and by 
Rottman^^ based on the Landau theory.^^ 

The new ingredient in the charge system, the long 
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range interaction, does not change the symmetry of the 
Landau free energy.^^'^^ StiU, the long range interaction 
can, in principle, change the order of the transition if 
the latter is dominated by fluctuations, making it devi- 
ate from the Landau theory scenario. While determining 
the specifics of the freezing transitions in this more gen- 
eral case is beyond the scope of the present work, one 
can make the following observations. 

(i) At the densities n = 1/3, 2/3, freezing occurs into 
one of the three degenerate x a/S configurations [see 
Fig. 6 (A) below]. In this state, electrons can occupy 
one out of the throe sublattices of the triangular lattice. 
We believe that this transition is of the first order. Our 
argument is based on the step-like singularity of the aver- 
age energy of the system at the transition observed dur- 
ing the MC dynamics, ^'^ and is corroborated by a fairly 
sharp step-like singularity in the MC conductivity (see 
below, Section IV), interpreted using the general connec- 
tion of the singularities in conductivity and in average 
energy.^'' Besides, for this transition a cubic invariant in 
the Landau free energy is allowed. ^"'^^'^"^ 

An alternative scenario for freezing at this density is 
a continuous transition of the q = 3 Potts universality 
class. ^^'^-'^ This possibility, in principle, cannot be ruled 
out based just on symmetry, since in two dimensions 
there are exceptions from the Landau theory, notably 
the g = 3,4 Potts models, for which the transition is 
second-order even though the cubic invariants exist. ^® We 
believe, however, that in the present case the existence 
of the cubic invariant triggers the first order transition. 
We also note that the experimental evidence for adsorb- 
tion of nitrogen molecules on graphite does not contra- 
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diet the first-order transition sccnario^^ even for fairly 
short-ranged quadrupolar molecular interactions. 

(u) At the density n = 1/2, the ground state is a (2 x 1) 
charge density wave of three possible orientations [see 
Fig. 6 (C)]. In the ordered state electrons occupy one out 
of the two sublattices, resulting in the six-fold degeneracy 
of the ground state. 

The nature of the phase transition for n = 1/2 is less 
clear. The classification^"'^^ suggests the g = 4 Potts uni- 
versality class with a second-order transition. However, 
the electron-hole symmetry at /it = line forbids the cu- 
bic invariant in the corresponding Landau free energy, 
putting this transition into the universality class of the 
Heisenberg model with cubic anisotropy.^^ The transi- 
tion in the latter model is still poorly imderstood,^® with 
both the continuous and the fluctuation-induced first or- 
der transition on the table. Recent study^* of the zero- 
field AIAFMwith finite-range interactions (beyond the 
nearest-neighbor) generally favors the first order tran- 
sition, while leaving room for an additional continuous 
transition at special values of couplings. At ^ 0, in the 
absence of electron-hole symmetry, the presence of the 
cubic invariant makes the first order transition scenario 
even more likely. 

Our numerical accuracy does not allow us to make a 
definite prediction. The observed kink in the conductiv- 
ity (Section IV below) is consistent with either the first 
or the second order transition, as is the singularity in 
the average energy. In general, freezing transition into 
the n = 1/2 ground state is less pronounced than that 
for n = 1/3, 2/3 states, since it is determined by the 
next-to-nearest neighbor coupling. 

(iii) For a generic density n we observe that at decreas- 
ing T the MC dynamics slows down, and all the charges 
eventually become immobile. The ground state in general 
looks disordered, and depends somewhat on the cooling 
history. The system configuration space appears to have 
many nearly degenerate minima, which complicates find- 
ing the true ground state numerically. While a disordered 
ground state for a continuous range of densities cannot be 
ruled out, we anticipate freezing into commensurate "epi- 
taxial solids" for at least some rational densities n = p/q 
with higher denominators, such as the striped ground 
state for the density n = 3/7 shown in Fig. 6 (D). These 
transitions would take place at ever smaller temperatures 
since they are governed by the couplings V{r) beyond 
nearest and next-to- nearest. For instance, for our model 
interaction (2) with d = 2a, freezing into the n = 3/7 
state occurs at T ~ 0.01V{a), and requires about 10 
hours of CPU time. We comment on possible freezing 
scenarios in Section X below. 



FIG. 3; Schematic phase diagram on the (/i, T) plane for a 
generic value of the screening length parameter d > 0. The 
dashed line marks the crossover between the correlated fluid 
and the disordered phase. The low temperature phases de- 
noted by a question mark could be either commensurate or 
disordered. 



IV. CHARGE TRANSPORT: ELECTRICAL 
CONDUCTIVITY 

Our interest in the hopping transport is two-fold. First, 
the conductivity is experimentally accessible in the dot 
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FIG. 4: (Color online) Temperature dependence of the zero 
bias dc conductivity cr{T). Shown are the curves a/n{l — n) 
for n = 1/3, n = 1/2, and for a typical intermediate density 
(taken here to be n = 3/7), obtained for 18 x 18 system with 
the screening length d = 2a. The asymptotic large tempera- 
ture behavior (21) is indicated by dashed line. The faint solid 
line, describing cr(r) for n = 1/2 on 12 x 12 patch, coincides 
with that for the 18 x 18 system, indicating the smallness of 
finite size effects. 



arrays."'^'''^^ Second, as we shall see below, the dc conduc- 
tivity is sensitive to the thermodynamic state, and its 
temperature and density dependence can thus be used to 
distinguish between different phases of the system. 

Below we use the MC procedure to compute the hop- 
ping conductivity in the presence of a small external elec- 
tric field. The MC conductivity ctmc and the one mea- 
sured in a real system are related as follows. Electron 
hopping between neighboring dots is assisted by some 
energy relaxation mechanism, such as phonons. The 
latter adds a temperature-dependent prefactor f(T) to 
the hopping rate, so that atotai = /(r)o-Mc(r). [In a 
simple model involving coupling to acoustic phonons, a 
Golden Rule calculation gives a power law /(T) cx T".] 
In this work, for simplicity, we shall ignore the system- 
dependent prefactor f{T), and focus on the MC conduc- 
tivity CTA/C- 

The temperature dependence of cta/c is shown in Fig. 4 
for several densities. The simulation was performed on 
a 18 x 18 patch using the charge-conserving dynamics 
(type B). The external field E, applied along the patch 
side, is chosen to be large enough to induce a current 
measurable in the presence of thermal fluctuations, and 
yet sufficiently small to ensure the linear response 

J-tE. (15) 

The data was obtained using E in the range 

Eac^i (10-2-5 X 10-2) y(a). 



We observed that the linearity holds for the field values 
much smaller than both the temperature and the next- 
to-nearest neighbor interaction, Ea <^ min{V^„„,r} 
[Eq. (18)]. 

The conductivity a, Eq. (15), is obtained using the 
MC dynamics of type B as follows. Since the field contri- 
bution to the potential difference between the adjacent 
sites separated by a = r^, is — Ea, we can incorporate 
the effect of the field by adding the quantity 



(ft - 9i)Ea 



(16) 



to the energy difference SElj in Eq. (8). 

Our simulation was carried out using the spin repre- 
sentation, as described in Sec. II C. In the spin language, 
using the Hamiltonian (11) and the temperature (14), 
10^ the term corresponding to Eq. (16) is (s^ — Sj)Ea. The 
data presented below was obtained using the field E ori- 
ented along the height of an elementary lattice triangle, 
so that |Ea| = Ea cos ^. The 'spin' current density was 
calculated as 



j = |(5s|acos(7r/6) 



(17) 



where J\f± is the number of hops along (against) the direc- 
tion of E during the MC run time, Af is the total number 
of MC trials at each temperature step, and Ss = ±2 is 
the spin change for each hop. The corresponding charge 
current is then jei = ^J- 

The MC conductivity as a function of temperature is 
displayed in Fig. 4 for several values of charge density. 
From the dependence cr(T) one can identify three tem- 
perature intervals with different behavior, corresponding 
to the high, low, and intermediate temperatures. Rele- 
vant temperature scales are approximately given by the 
nearest neighbor and next-to-nearest neighbor interac- 
tion strength: 



(18) 



In our MC simulation, for d = 2a, the value of Vnnn, 
given by the interaction across the main diagonal of a 
rhombus, was about 0.3y„„. 

The simplest to understand is the high temperature 
behavior T > V^„, corresponding to a disordered phase 
in which conductivity takes place via uncorrelated hops 
of individual electrons. Conductivity in this phase 
has a simple temperature and density dependence [see 
Eq. (21)]. 

In the opposite limit, at T ^ 0, the system freezes 
into the ground state configuration (solid phase), and 
the conductivity vanishes. The freezing transition is en- 
tirely due to the longer range coupling, such as Vnnm 
since for purely nearest neighbor coupling, realized in 
the AIAFM model, the system does not exhibit a phase 
transition and is characterized by finite entropy even at 
T = 0. Thus the upper temperature scale for the sohd 
phase can be estimated as T ^ Vmm- The ground state 
depends on the density n in a complicated way. Near 
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rational n the ground state is commensurate, while at 
a generic n the state is probably incommensurate. The 
freezing temperature Tc is also a function of the density 
n. As the temperature approaches the freezing temper- 
ature Tc, the conductivity vanishes (Fig. 4). While the 
behavior <j{T) near T — Tc appears to be singular (see 
Fig. 4) , we were not able to extract the exact form of this 
singularity directly from our simulation. However, some 
information about the singularity in conductivity can be 
obtained from general arguments^^ relating it with the 
distribution of energies studied in Sec. III. We also ob- 
serve singularities near rational n in the conductivity de- 
pendence on the density, a{n). They are indicated by 
arrows in Fig. 5. 

Finally, at densities 1/3 < n < 2/3 there is an inter- 
esting intermediate temperature interval 

Vnnn <T<Vnn (19) 

in which the conductivity is finite and reaches maximum 
as a function of temperature (Fig. 4). Transport at these 
temperatures is of a collective character [correlated fluid 
phase) owing to strong short-range correlations between 
charges which constrain individual charge movements. 
Due to the frustration, the onset of short-range ordering 
docs not lead to charges freezing. The system appears to 
possess a sufficient amount of residual entropy to allow 
for finite conductivity. 

An independent consistency check for the MC dynam- 
ics is provided by the fluctuation-dissipation theorem, 
relating current fluctuations to conductivity: 



(20) 



Our simulations is found to be in accord with (20) in all 
temperature regions, which ensures that the MC conduc- 
tivity indeed describes transport in the linear response 
regime. (Some deviations from (20) were observed very 
close to the freezing point, where a becomes very small.) 
In the correlated phase we explicitly evaluate the inte- 
gral in the left hand side of (20) by numerically averag- 
ing the product j^(t -f T)jy{t) over r and t. The result is 
compared with the conductivity obtained directly from 
Eqs. (15), (17) to make sure that during an MC run the 
system has enough time to reach equilibrium. 

At large temperature T 3> V{a), one can evaluate the 
left hand side of the fluctuation-dissipation relation (20) 
explicitly and find the universal high temperature asymp- 
totic behavior of the conductivity, 

2 "(1 - 



a = a 



T 



(21) 



To obtain (21) we note that for a high enough tempera- 
ture the current is delta-correlated in time, 



(22) 



The mean square (j^) can be obtained by summing the 
probabilities of possible MC moves: 

(f) = ^ • 2n(l - n) ■ {6sf ■ (a cos J) ' • (23) 



The factor | comes about because in the field E aligned 
along the height of the elementary lattice triangle, only 
four out of possible six bond directions contribute to con- 
ductivity. The second factor, 2n(l — n), describes the 
probability to select two adjacent sites occupied by an 
electron and a hole. The change of occupancy per MC 
move is Ss — ±2, since Si = ±1. The expression (23) does 
not depend on temperature since all the hops are equally 
probable and uncorrelated, Wi^j = Wi^i = w = 1/2 
at r ^ V{a). The conductivity tensor is isotropic, 
= f(5^i^, as can be explicitly checked by doing a 
similar high temperature calculation for an arbitrary ori- 
entation of the field E. Fig. 4 shows that the results of 
our MC dynamics are consistent with the high tempera- 
ture behavior (21). 

The change in conductivity behavior while cooling 
down from the disordered into the correlated fluid and 
solid phases can also be seen in Fig. 5. Here we plot 
the conductivity cr, scaled by a? /T, as a function of the 
electron density n for several temperatures. The results, 
shown for the densities < n < 1/2, can be extended 
to n > 1/2 using the electron-hole symmetry n <-> 1 — n. 
The high temperature curve is clearly consistent with the 
n{\—n) dependence (21). The low temperature curves in- 
dicate that the conductivity vanishes more quickly near 
the densities of a simple fraction form (n = 1/4, 1/3, 
1/2). This corresponds to freezing of the system into a 
commensurate state at these values of n. The commen- 
surate states at n = 1/3 and n = 1/2 are shown in Fig. 6 
(panels A and C respectively). 
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FIG. 5: (Color online) The conductivity a, scaled by a? /T, is 
shown as a function of electron density for several tempera- 
tures. The temperature values are given in the units of V[a)\ 
the interaction V{r) of the form (2) was used with the screen- 
ing length d = 2a. Arrows mark the features corresponding 
to the freezing phase transitions at n — 1/4, 1/3, 1/2. Dashed 
line corresponds to the high temperature limit described by 
Eq. (21). 
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FIG. 6: (Color online) A: Commensurate \/3 x ^/3 charge 
ground state at the density n — 1/3. Electrons are are rep- 
resented by small dots on the triangular array sites. Pairing 
of the triangles is revealed by erasing all frustrated bonds 
connecting the sites with equal occupancy, as described in 
Section V. B: Typical charge configuration for n — 1/3 + e. 
There are three excess charges in the system (e = 3/144) hop- 
ping over the honeycomb network of unoccupied sites in the 
state shown in panel A. The excess carriers, dilute at e ^ 1, 
are moving nearly independently on the frozen ^/S x y/S state 
background. C: Ground state for n = 1/2. Shown are the 
two charge density wave domains characterized by different 
slope orientations [Eq. (26)]. D: Example of a ground state 
for a commensurate density with a higher denominator, here 
n — 3/7. E: Typical charge configuration obtained in a sim- 
ulation for n = 1/2 in the correlated fluid phase. Two el- 
ementary '3d cubic' cells corresponding to a free charge are 
marked. F: Typical charge configuration for n = 1/2 at tem- 
perature somewhat higher than in panel C. In this case, there 
are several topological defects (unpaired triangles) present in 
the system. These defects can be interpreted as dislocations 
of the height field (see Sec. V) . 



At densities near these simple fractional values, the 
system conducts via hops of excess electrons or holes 
moving in the frozen crystalline background. Such a sit- 
uation is depicted in Fig. 6 (panel B) for the case of 
n = 1/3 -I- e. Since the conductivity is proportional to 
the excess charge density e, we expect a{n) to have cusps 
near simple fractional densities. Such cusps are indeed 
seen in Fig. 5 near n = 1/4, 1/3, 1/2. 



V. THE HEIGHT VARIABLE 

A. The ground state 

Here we attempt to understand the ordering in the 
ground state and in the correlated fluid by employing the 
notion of the height field order parameter, originally in- 
troduced in the context of the AIAFM problem by Blote 
et a/.^^'^* In our model (1), (2), the latter problem cor- 
responds to the limit d <C a dominated by the near- 
est neighbor interactions. The essential features of the 
ground state can be understood by considering energy 
minimization for individual plaquettes (triangles). With 
only two kinds of charges it is impossible to avoid the 
frustrated bonds between the like charges, with at least 
one such bond per triangle. In the ground state, the num- 
ber of such bonds must be as small as possible. That can 
be achieved by pairing the neighboring triangles in such 
a way that each pair shares a frustrated bond. One can 
see that the charge state corresponding to all triangles 
paired, while providing the absolute energy minimum, is 
not unique. On the contrary, the pairing condition leaves 
plenty of freedom in the charge configuration, character- 
ized by extensive entropy (finite entropy-to-area ratio). 

The configurations permissible by the pairing condi- 
tion have a simple geometric interpetation. After erasing 
all frustrated bonds, one obtains a rhombic tiling of the 
triangular lattice, ^^'^^ as shown in Fig. 6 (panels A - E). 
It is convenient to view such a tiling as a 2d surface in a 
3d cubic crystal projected along the (111) axis on a per- 
pendicular plane. After undoing the projection by lifting 
the 2d configuration of rhombi in the 3d space, each site 
Ti acquires a scalar variable h{ri) which is the height of 
the lifted surface in a 3d cubic crystal. The field h takes 
values which are multiples of the distance between the 
cubic crystal planes, 

where £ is the main diagonal of the unit cell in the cubic 
crystal. The mapping onto a continuous height surface 
exists only for the densities 1/3 < n < 2/3, which we 
focus on hereafter, since for n outside this interval the 
tiling contains voids. 

As we discussed above, the long range character of the 
interaction (2) lifts the ground state degeneracy, leading 
to freezing into specific ground states. To use the height 
variable to characterize the ordering, we define the slope 
t of the height surface. For that we consider the surface 
normal vector 

3 

m = A^^Yl ' A = ri+r2+r3 (25) 

where ri, r2, are the numbers of rhombi formed by 
(62,63), (63,61), (61,62), which add up to the total area 



9 



A = N"^ ■ The slope t is given by the projection 

t = m-(z-m)z, z = (ei + 62 + 63) . (26) 

Here the vectors e^, s = 1,2,3 form the basis of the 
auxiliary 3d crystal (Fig. 10 in the Appendix), and 
give the numbers of 3d crystal faces normal to the vectors 
e^. Eq. (26) defines an average slope of the configuration. 
For a 'smooth' surface, the vector t(r) is proportional to 
the local height gradient, t(r) cx {dxh,dyh). 

When the interactions are of the nearest neighbor kind 
(d <C a), the height surface fiuctuates freely even at 
small T. We observe a similar behavior at finite tem- 
peratures (above freezing) for the system with the long 
range interaction (2). The numbers r^, s = 1, 2, 3, for an 
N X N patch are equal on average, fiuctuating around 
the mean value of -/V^/3, and thus the average slope is 
t = 0. As T ^ 0, for the studied densities n = 1/3, 
2/3, and n = 1/2, the interaction (2) leads to freezing 
into commensurate ground states characterized by spe- 
cific discrete slopes and a non-extensive entropy. 

B. The correlated fluid 

The height field h is defined globally and uniquely (mod- 
ulo an additive constant and an overall sign) for any of 
the degenerate ground states of the AIAFM model. For 
T > 0, however, due to thermal fluctuations, some of the 
triangles are left unpaired. Such unpaired triangles rep- 
resent topological defects of the height field. This can be 
seen most clearly at low T, when the defects are dilute, 
and the height field can be constructed locally around 
each defect. It then turns out that the height field, con- 
sidered on a closed loop surrounding a defect, is not a 
single valued function. 

The topological charge assignment in this situation is 
facilitated by interpreting the defects as screw disloca- 
tions, centered at the unpaired triangles. Positions of the 
dislocations are seen as large (2 x 2) triangles in Fig. 6 
(F). After integrating V/i over a contour enclosing a sin- 
gle triangle, one obtains a positive or negative mismatch 
in the height h equal to the dislocation Burgers vector 

h^= <j> d^hdx^^±2l^±6b. (27) 

The topological charge algebra is indeed identical to that 
of dislocations [Z). Following a loop around two dis- 
locations of opposite sign (27), gives net zero charge, 
§ dihdxi = 0. The topological nature of the defects con- 
strains the dynamics. During the MC simulation, the 
dislocations originate and disappear only in pairs. 

The minimal energy cost required to create a defect 
can be estimated as the energy for unparing of two tri- 
angles. Each unpaired triangle adds one more frustrated 
bond which can be shared by neighboring triangles [see 
Fig. 6 (F)]. The corresponding energy cost is 2V{a) per 



unpaired triangle. At sufficiently low temperatures the 
probability of creating a pair of defects is therefore ex- 
ponentially suppressed, P ~ e-4y(a)/T_ ^j^^ other 
hand, at a high temperature T ^ V{a) the defects are so 
abundant that the height field ceases to exist even as a 
local notion. In this high temperature state, marked dis- 
ordered phase in Fig. 3, electron hopping is uncorrelated, 
with the interaction enforcing single or zero occupancy, 
and otherwise playing no role. 

At lower temperatures, T < V{a), with the fugacity of 
a single defect decreasing at least as e~^^^°)/-^, the de- 
fects quickly become dilute. In this case, the height field 
is defined locally in the entire plane, except the vicinity 
of the defects. The collective behavior of defects in the 
correlated phase, such as the Kosterlitz-Thouless unbind- 
ing transition, is controlled by entropic effects which will 
be analyzed below in Section VII. 

The charge rearrangements that do not produce or de- 
stroy defects cost no energy in the AIAFM model limit 
d a. In the case of a long range interaction, such re- 
arrangements generally cost finite energy which is deter- 
mined by the Vnnn interaction strength. An example of a 
typical single electron move that preserves the structure 
of rhombic tiling is shown in Fig. 7. We find that the state 
with a height field, despite being constrained by the tri- 
angle pairing condition, allows for sufficiently large num- 
ber of movable charges which can lead to macroscopic 
charge rearrangements which do not produce defect pairs. 
The MC conductivity appears to be related to these re- 
arrangements, hereafter referred to as free charges. As 
Fig. 7 illustrates, the moves associated with free charges 
can be interpreted geometrically as changing the 3c?-lifted 
tiling by adding or removing two adjacent cubes. This 
operation does not introduce defects or discontinuities in 
the lifted surface, preserving its 3d continuity. 

In our MC study, with electron densities 1/3 < n < 
2/3, we found that there is a temperature interval in 
which the number of topological defects (unpaired tri- 
angles) is small, and thus the height order parameter is 
well defined. Simultaneously, the density of the free 
charges (see Fig. 7) was observed to be large. During MC 
simulations we accumulate a histogram for n/ and find 
that this is a nonconserved quantity with a broad Gaus- 




FIG. 7: (Color online) The concept of free charges in the cor- 
related phase, which can move without changing the number 
of frustrated bonds, is illustrated. Above: An example of a 
free charge move; Below: The corresponding change in the 
rhombic tiling. 
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sian distribution. The mean n f is of the order of 10—50 % 
of the total electron density n for the temperatures 

Vnnn < T < V^n (28) 

between the nearest neighbor interaction V^„ which en- 
forces pairing correlations, and the next-nearest interac- 
tion Vnnn which controls freezing at low temperature. 

The correlated nature of the charge state at these tem- 
peratures, owing to the continuous height variable con- 
straint, leads to a peculiar picture of charge flow. The 
notion of individual electron hopping has to be replaced 
by a more adequate picture involving fluctuations of the 
height field giving rise to charge movement. We observe 
that the Ohmic conductivity remains finite (Fig. 4) in 
the correlated fluid temperature interval, while the topo- 
logical defects freeze out. It appears that the height fluc- 
tuations on their own are sufficient for charge transport 
and conductivity. The elementary height fluctuations, 
embodied in the notion of free charges (Fig. 7), represent 
MC moves of an effective nonconserving (type A) dy- 
namics of the 2d surface. Based on this idea, in Section 
VIII we shall develop a continual approach to describe 
the dynamics in the correlated state. 

One may expect that the free charges, although nom- 
inally allowed to only move back and forth (Fig. 7), can 
propagate over the whole system. Microscopically this 
happens since a move of one free charge unlocks subse- 
quent moves of other free charges. This picture is con- 
sistent with our observations based on MC dynamics: In 
Sec. VI we find a non-vanishing contribution to the con- 
ductivity in the absence of dislocations. The dislocation- 
less conductivity points to the importance of free charges 
in transport. However, more work will be needed to ad- 
dress other properties of free charges dynamics, such as 
the microscopic transport mechanism and crgodicity. 

To assess the relevance of this conductivity mecha- 
nism for real systems one needs to understand several 
issues, the most important one being the role of disor- 
der. Although our MC study of the disordered problem 
was not too extensive and, in partricular, restricted to 
not too low temperatures, it allows to draw some assur- 
ing conclusions. In general, for weak disorder we observe 
no change in the qualitative features of the dynamics. 
We studied the dc conductivity in presence of a finite 
amount of disorder, modeled by the random potential 
(/>(ri) in the Hamiltonian (11). The MC simulation used 
statistically uncorrelated random (p{vi) with a uniform 
distribution in the interval — 0o < < 0o, where 

(/•o ~ Kinra = V(\/3a). The only significant difference 
observed was a relatively more slow time averaging in the 
presence of disorder, leading to enhanced fluctuations of 
the conductivity recorded as a function of temperature. 
Otherwise, the qualitative behavior of the conductivity 
was found to be the same as in the clean system. In par- 
ticular, the zero bias conductivity remains finite in the 
correlated fluid phase in the presence of disorder. We at- 
tribute the robustness of conductivity against weak disor- 
der to the nonconserving dynamics of the height variable, 



which is not pinned by the disorder in the temperature 
interval of interest. 



VI. DISLOCATION-MEDIATED VERSUS 
HEIGHT FIELD-MEDIATED CONDUCTIVITY 

Here we investigate in detail the effect of dislocations 
on electron transport in the correlated phase. As 
noted above, certain charge moves, associated with "free 
charges" (Fig. 7) have relatively low energy cost. These 
moves correspond to height field fluctuations which do 
not produce topological defects. The presence of the pair- 
ing correlations described by the height variable and the 
dynamics of free charges associated with these correla- 
tions, poses an interesting question regarding the relative 
contribution of the free charges and dislocations to the 
transport. 

One approach to understand the role of dislocations 
would be to forbid entirely the charge hops that introduce 
new dislocation pairs and study MC conductivity in such 
a system. We have tried to modify the MC simulation in 
this way, and found that the dislocationless system ex- 
hibits finite conductivity solely due to free charges, indi- 
cating that the dislocations are not essential for conduc- 
tivity. This method, however, docs not allow to compare 
the role of dislocations and free charges quantitatively, 
since the rule introduced to eliminate dislocations alters 
the dynamics in an unc;ontrollablc way. In fact, since the 
definition of conductivity changes in the dislocationless 
MC, although we indeed observe a nonzero conductivity, 
it is difficult to compare the resultant "dislocationless" 
conductivity with that of the original problem. Instead, 
we adopt a different, more gentle approach. We employ 
the same dynamics as above (conserving, or type B), in 
which we detect dislocations and analyze partial conduc- 
tivities ap in the presence of p = 0, 1, 2, .. dislocation pairs 
(Fig. 8). 

Let us outline the corresponding MC procedure. We 
introduce a small bias E across the finite system with 
periodic boundary conditions, and define the partial con- 
ductivities (jp similarly to the total conductivity calcu- 
lated as described in Section IV, Eq. (15): 

jp = <Tp^- (29) 

The partial currents jp at a fixed number p of the dis- 
location pairs are defined similarly to the total c;urrcnt, 
Eq. (17), with the numbers of hops N±, N ^ U+ + 
replaced by the p-dependent Mp^ and Np = Npj^ +Np_. 
Here Mp is the total number of MC trials, or attempted 
hops that may or may not lead to real hops, in which a 
pair of randomly selected neighboring sites is such that 
the charge hop between them would preserve the number 
j3 of dislocation pairs. Accordingly, A/^^ andA/^_ descibe 
charge hop trials along and opposite to the applied field 
with a fixed number p of dislocations. The actual hop 
is then attempted with a Boltzmann probability (7), as 
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FIG. 8; (Color online) Partial conductivities o-p, p = 0, 5, 
which account for MC moves in the presence of p dislocation 
pairs (see text), for the 12x 12 patch, density n = i, screening 
length d = 2a, Eq. (2). The dislocationless conductivity oq 
approaches zero continuously at the freezing transition, while 
the conductivities a^, p ^ Q remain finite. Separately shown 
are the conductivities cr| and ay that correspond to the MC 
steps changing the number of dislocation pairs. (Note the 
negative sign ay < 0.) 




Temperature, T 

FIG. 9: (Color online) Partial attempt numbers scaled by 
N define normalized weights Np/N , p =0,...,5, ti ii for the 
partial conductivities of Fig. 8. The total N [Eq. (30)] was 
approximately constant, varying by less than 5% in the whole 
temperature range. With the A/" « 3 • 10* per temperature 
step, the whole simulation took about 15 hours on a 1.25 GHz 
PowerPC. 

discussed in Section II. In our simulation we studied the 
system 12 x 12 at density n = i. 

The partial conductivities Up, p = 0, 1, 2, describing 
the contribution to transport of charges hopping in the 
presence of the 2p dislocations, need to be complemented 
by two additional parts, (T| and (T|, describing the con- 
duction processes which alter the number of dislocations. 
Indeed, a significant fraction of the occupied-empty site 
pairs chosen in the MC simulation are such that upon a 
charge hop between them the dislocation number would 



either increase or decrease. Such hops, if selected by the 
Boltzmann probabilities, can also carry current in the 
presence of the applied bias (as illustrated in Fig. 8). 
With the corresponding numbers of the MC trials de- 
noted by A/| and M[ , the total number of attempts is 

U = No+Ni+N2 + (30) 

The partial MC attempt numbers Mp, A/|, A/| tempera- 
ture dependence is summarized in Fig. 9. 

As Figs. 8 and 9 illustrate, the average number of dis- 
locations in the 12 x 12 patch evolves from about 5 at 
T w Q.lV{a) to zero in the vicinity of the freezing tran- 
sition. As the temperature lowers, most of the time the 
system attempts to create dislocation pairs (A^ is dom- 
inated by A/j), however, almost all of these attempts 
are discarded due to their exponentially low Boltzmann 
weight, resulting in the negligible current j| and partial 
conductivity cr| (see Fig. 8). 

We note a drastic difference between the dislocation- 
less conductivity uo and the partial conductivities ap^Q 
near the freezing transition at Tc ~ Q.2V{a). Whereas 
the former approaches zero continuously as T —> Tc, the 
latter have an apparent step-like discontinuity. However, 
the contributions of partial conductivities Up, p 7^ 0, to 
the total conductivity 

(j = N-^ MpUp (31) 

p=0,l,2,...,T,i 

are small, since the weights Np^o/N drop very quickly 
as T approaches Tc (Fig. 9). In a relatively small 12 x 12 
system studied here, the dislocations are almost always 
totally absent near Tc. In this case the total conductivity 
(31) is dominated by ctq (Fig. 8), due to small Np=^o/N. 
This is consistent with the observation that the total con- 
ductivity decreases to zero in a continuous fashion near 
Tc, similar to ctq, with the step-like contribution of Up^Q 
being inessential. 

One cannot exclude, however, a different regime in 
a large system, where the discontinuous part Up^Q of 
the total conductivity (31) can become dominant in the 
thermodynamic limit, due to dislocation number growing 
with system size. Even in this case, we expect the dislo- 
cationless contribution to the conductivity to remain sig- 
nificant. In that regard, we note an approximately con- 
stant increment in conductivity Apcr = ap+i — Gp ~ const 
when the number of dislocations increases by one, p — > 
p + 1 , observed at temperatures not too close to freezing 
(Fig. 8). The approximately p- independent ApCr provides 
an estimate of the conductivity per dislocation, suggest- 
ing that in the dilute regime the dislocations contribute 
to transport approximately independently. However, the 
dislocationless part ctq is a few times larger than Aptr, 
providing a constant offset to the Op vs. p dependence. 
It is thus not unconceivable that, if the dislocations are 
suppressed, e.g. due to Kosterlitz-Thouless transition, 
or otherwise, the conductivity would remain finite, dom- 
inated by the dislocationless contribution. 
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This conclusion is reinforced by the resuhs of MC sim- 
ulation with the rules modified so that the dislocations 
are totally excluded, as discussed above. In this case, we 
find finite conductvity, behaving as a function of temper- 
ture at T > Tc in a way similar to ao- Interestingly, the 
freezing transition is unaffected by this alteration, with 
the same value of Tc observed under changed MC rules. 
These observations, in our view, leave no doubt that the 
conductivity docs not entirely depend on dislocations. 
The transport of free charges, responsible for the dislo- 
cationless contribution ctq, gives rise to finite conductivity 
independent of dslocations. 

Interestingly, the partial conductivity due to annihila- 
tion of the topological defect pairs is negative, ai < 0. 
While we do not have a complete understanding of this 
observation, we mention one possible explanation which 
involves a nonlinear effect. We suppose that a large 
enough electric field can perturb the system so that its 
relaxation back to equilibrium will be accompanied by 
release of a charge in the opposite direction. However, 
our MC study indicates that the relative weight of such 
processes, and thus its contribution to the conductivity, 
is small near the freezing transition (Fig. 9). Therefore, 
in our view, the associated nonlinear conductivity is not 
essential for the interpretation of the numerical data. 

The described analysis of partial conductivities does 
not reveal the microscopic mechanism by which disloca- 
tions facilitate transport. One possibility is that the con- 
duction is due to the motion of dislocations themselves 
(since they are charged objects). Another possibility is 
that the mere presence of dislocations facilitates the hops 
of free charges in their vicinity and results in a finite 
conductivity. While our MC study does not provide a 
definite answer, it is possible that the two alternatives 
are not unrelated, since the motion of a dislocation can 
happen due to a free charge movement right next to it. 

To summarize, we have shown that the dislocation 
pairs facilitate conductivity, and have compared their 
contribution with dislocationless conductivity due to free 
charges. While in the model analyzed here there is no 
obvious way to completely separate these contributions, 
due to the absence of a dislocation binding transition, we 
conclude that both contributions are present, with their 
relative importance depending on the dislocation density. 



VII. FLUCTUATIONS IN THE CORRELATED 
PHASE: DISLOCATIONS UNBINDING 

The height field describes the effect of frustration on 
charge dynamics by providing a nonlocal change of vari- 
ables that helps to keep track of the local correlations 
between charges. However, as we found above, in gen- 
eral the height variable cannot be globally defined due 
to the presence of dislocations. Nonetheless, at relatively 
low temperature, when dislocations are dilute, the height 
field provides a useful description on a local level. 

Here we adopt the view that the dynamics of the height 



field is simpler than the underlying charge dynamics. In- 
deed, while the latter is strongly constrained, the fluc- 
tuations of the height field appear to be Gaussian. The 
height variable also provides a natural continual descrip- 
tion employed to study dislocation unbinding (this Sec- 
tion), roughening transition (Section VIII), the dynamics 
(Section IX) as well as the freezing into commensurate 
states (Sec. IIIB above). 

The MC dynamics in the correlated state can be inter- 
preted as the height field fluctuations. These fluctuations 
dominate in the temperature interval (28) where the dis- 
locations are dilute. Thus here we consider the height 
field ignoring for some time the effect of dislocations. We 
study the height fluctuations numerically and find that 
in the continuum limit r ^ a they are described by the 
partition function of the form 

Zo = ym(r) e" (32) 

where k is an effective stiffness of the height surface, with 
temperature incorporated into k. 

The Gaussian partition function (32) has been intro- 
duced by Blote and Hilhorst,^'' and by Niehnuis et al.^* 
for the height field in the AIAFM model at T ^ 0. The 
corresponding stiffness has been obtained from the exact 
solution^'5'1^: 

7r 

kaiafm = • (33) 

[Note that, instead of the stiffness k, Ref. 18 employs the 
"renormalized temperature" Tr = 2-jT/Kb^, with Tr = 18 
for the AIAFM .] The equilibrium height fluctuations 
have been studied numerically''''''^^ in the case of the near- 
est neighbor interactions, confirming the result (33) for 
the AIAFM . 

Here we extend these results to the model with long 
range interaction. We study the fluctuations of the height 
field numerically, by assigning the heights h{ri) to the 
triangular lattice sites according to the procedure de- 
scribed in Appendix. The results of this study can be 
summarized as follows. 

(i) At fixed temperature, we accumulate the his- 
tograms 'P[h{ri)] of height values for a set of points {vi} 
in the N x N patch. We find that the fluctuations are 
Gaussian, 

\ogP[h{ri)] oc -h\ri) , (34) 

independently for each point r^. This suggests that a 
partition fimction of the form (32) can be employed. 

(ii) To calculate the effective stiffness for a system 
with the long-ranged interaction (2), we fit the two-point 
height correlator to the one following from Eq. (32): 

{ihir,)-h{r,)r)^l-lJ-^, (35) 

with the length vq of the order of the lattice constant a. 
For this analysis we generate random MC configurations. 
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selecting the height surfaces of zero average slope, defined 
as in Eq. (26) above. 

We tested this procedure by calculating the stiffness in 
the nearest neighbor interaction limit (i<Caatn = l/2. 
The MC dynamics recovers the value kaiafm, Eq. (33), 
within a few percent accuracy. With the long range in- 
teraction, the relationship (35) still holds, allowing us 
to determine k. We systematically find the values of k 
below KAiAFM, Eq. (33). 

(iii) The measured stiffness k can be used to assess 
the possibility of a dislocation-binding transition of the 
Kosterlitz-Thouless kind.''^''^'' If exist, such a transition 
between the disordered and correlated fluid phases would 
realize34.35 the KTHNY 2d melting scenario. Dislocation 
binding takes place above the universal threshold value 

K > kkt = To- , (36) 

"a 

with &A the Burgers vector of the dislocation. Using the 
value from Eq. (27), one obtains kkt = 2k;aiafm- Thus 
the KT transition is absent in the zero-field AIAFM }^ 

We find that this conclusion remains valid for the 
long range interaction. For the interaction (2), the stiff- 
ness decreases as a function of d. For typical d used in 
this work, stffness values are a few times smaller than 
AIAFM • Although such a behavior may seem counterin- 
tuitive (one could expect the system to become stiffer 
for an interaction of longer range), it is in fact con- 
sistent with the scaling arguments presented below in 
Sec. VIII B. It is also compatible with the observation^* 
that the ferromagnetic second-neighbor coupling causes 
increase of stiffness. In agreement with this, adding 
longer-range antiferromagnetic couplings should lead to 
the decrease of k. 

We conclude that the dislocations are always unbound 
in our system. This is indicated by a crossover (rather 
than a sharp transition) between the disordered and cor- 
related phases in the phase diagram, Fig. 3. Thus, in the 
thermodynamic limit, the dislocation pairs are present 
at any temperature. While their concentration may be 
small, it is nonzero, and, strictly speaking, the height 
field is not well-defined in an infinite system. However, 
in a finite system we often observe that dislocations are 
absent at low temperature, which justifies the height field 
as a local notion. 



VIII. THE ROLE OF HEIGHT FLUCTUATIONS 
NEAR FREEZING TRANSITION 

Although the height surface fiuctuates freely in the cor- 
related phase, the stiffness values k obtained from the 
MC dynamics are found to be below the bound^* neces- 
sary for the fluctuation-induced roughening transition.''® 
In this respect, the case of the long-range interaction 
(2) appears to be qualitatively similar to that of the 
AIAFM in the small fleld.^*'^'* The purpose of the present 
Section is to rationalize this similarity, and quantify 



the reduction of stiffness (as compared to that of the 
AIAFM ) by employing renormalization group methods. 

Below we construct the free energy for the charge den- 
stity and the height field, with a coupling of a sine- 
Gordon form, and present its scaling analysis. We ob- 
tain a bound on the effective stiffness k, Eq. (54), and, 
by comparing it to that found from the MC dynamics, 
determine that the continuous roughening transition does 
not take place. Rather, as discussed in Section III, freez- 
ing occurs via a finite order transition. The bound (54) 
indicates that the height fluctuations are irrelevant, and, 
as a result, allows us to rule out both instances of the con- 
tinuous transitions (dislocation-unbinding, and roughen- 
ing). Subsequently, in Section IX, we use the developed 
formalism to study equilibrium current fluctuations and 
conductivity. 

A. Free energy 

The correlated state is described in terms of the height 
variable h{r) and charge density n(r). Microscopically 
the state of the system is defined uniquely by specifying 
either the former or the latter. However, here we demon- 
strate by scaling analysis that at large length scales the 
height and density fluctuations decouple. This obeserva- 
tion allows us to employ a partition function with h and 
n as independent variables, 

Zcorr = y"m(r)2?n(r) e-^[^'"l . (37) 

Here the free energy 

J^[h, n]=Th+Tn+ Tint (38) 

is a sum of the contributions of the height field, the charge 
density, and their interaction. As before, the temperature 

is incorporated into T . 

The phenomenological free energy J^h has the form 

Th = y"d'r(^|(Vr/i)2+5cos^+/cos^^ , (39) 

with K the stiffness. Here the term cos assigns higher 
statistical weight to the field configurations that pass 
through the points of the 3d cubic lattice, with the period 
h in height given by Eq. (24). The third term in (39) de- 
scribes the coupling of charge to the chemical potential, 
in our case realized as the gate voltage, 

/ = CneV^ , (40) 

with n the average charge density and C ~ T~^. The 
coupling of the form cos with the period 26 in height, 
arises because opposite charges occupy two different 
sublattices of the 3d lattice, alternating in the height 
direction. -'^^ The third term in Eq.(39) is always more 
relevant than the second one in the sense of scaling. 
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Consider now the charge density nr. Electrons inter- 
act with each other and with an external electrostatic 
potential $'"^*(r), which gives 

jc-„ = i jcfrcfr'nrUr-r'Ur' + jd^rnr<^'"'\r) (41) 

with the interaction Ur-r' of the form (2). 

The coupling between the density and the height vari- 
able should be periodic in 2b for the same reason as in 
Eqs. (39) and (40), with n replaced by the fluctuating 
density: 

^int ~ ^ j • ('^^) 
Scaling analysis of the problem (38) is presented below. 

B. Scaling analysis 

The nonlinear terms in Eqs. (39), (42) acquire nontrivial 
scaling dimensions due to fluctuations. At one loop, the 
anomalous dimensions of the couplings /, A and g are 



Xf = X\ = 



X„ 



2k* 

K 

Sk* 



where we define 



K = 



TT 

862 ■ 



(43) 

(44) 

(45) 



The term / cos ^ becomes relevant when Xf < 2, or 

K > Kf = K* . (46) 

The coupling g becomes relevant at a larger stiffness; 

K > Kg = Ak* . (47) 

The coupling A between the density and height fields 
is relevant when the scaling dimension of the gradient 
term is larger than the sum of scaling dimensions of A 
and n. The field n is not renormalized since its free en- 
ergy is quadratic. From (41), the bare dimension of n is 
Xn = 3/2. Hence the height-density interaction is rele- 
vant when 



X^ 



Xn < 2 , or K > K\ = 4k* 



(48) 



The nonlinear terms in (39) do not renormalize the 
stiffness k at one loop. However, the stiffness may obtain 
a one loop correction 5k dnc to the interaction with the 
density, Eq. (42). Below we consider such a possibility. 

Since the free energy ^„ is Gaussian, we integrate out 
the field Ur in (37) and obtain an effective free energy for 
the height field: 



2 ^ 

k 



(49) 



Consider the case of the screened Coulomb interaction of 
the form (2), 



27r 



(50) 



At length scales larger than the screening length, /cc? ^ 1, 
1 1 



C/k 27r 



{A Bk) 



(51) 



with A = 1/d, B = 1/2. For the unscreened in- 
teraction, (C/k)^^ is of the form identical to (51) with 
A = 0, B = 1. The local term A/2tt corrects the cou- 
pling constant g and does not contribute to the stiff- 
ness. The correction to stiffness 5k arises from the 
nonlocal part of (49) with U^^ Bk/2iT. Expanding 
cos-Kh/b = [e^-^^l^ + e-''^^/^) /2, we note that the terms 
^^iTrh/b-^^^^piTTh/b-^^ + c.c. havc a larger scaling dimension 
than the cross terms, (e-*''''/'')_k(e''''*/'')k + c.c. . Ne- 
glecting the former and expanding the latter, we arrive 
at the following stiffness correction: 



5^ = - ^ ^ JkK*k'^h-\Ji\^ 



(52) 



where we introduced a nonlocal coupling Jj, = 2i?A^/|k|. 
Note that the screening length d does not contribute to 
the stiffness correction, affecting only the value of B. The 
negative sign of the stiffness correction (52) is consistent 
with the observed stiffness reduction in MC simulation 
with long range interaction, compared to kaiafm- 

The coupling J is renormalized via integrating out the 
height fluctuations /ik with l/li < k < l/lo- 



t. In \ I 



1— 4k*/k 



(53) 



The stiffness correction (52) is relevant only when the 
condition (48) holds. Therefore, there is no qualitative 
effect on the Gaussian behavior (32). 

We point out that the negative contribution to stiff- 
ness, Eq. (52), can be understood as an effect of long 
range interaction as follows. In the 3d cubic crystal pic- 
ture, the charges of plus and minus sign reside on the 
even and odd sublattices. With respect to the (111) di- 
rection used to introduce the height variable, these sub- 
lattices define alternating stacks of planes normal to the 
(111) axis, each filled with charges of one sign. Now, ev- 
ery surface with a small gradient of the height field can 
be associated with a set of terraces of alternating total 
charge and of width inversely proportional to the gradi- 
ent of the height field. Since the long range Coulomb 
interaction favors the states with plus and minus charges 
well-mixed, one expects it to favor more narrow terraces 
compared to wider terraces. This means that the long 
range repulsive interaction contribution to the energetics 
of the height field must be negative of (Vh)"^, so that the 
favored states have the maximal possible height gradient. 
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This observation is consistent with both the negative sign 
of the result (52), and with the stripe- like character of 
the ground state at n = 1/2 [Fig. 6(C), single domain], 
stabilized by the long-range repulsion. 



C. Bounds on stiffness 

To summarize the above results, the scaling relations set 
bounds on the effective stiffness k in the freely fluctuating 
height model, Eq. (32). The Gaussian behavior is stable 
when the nonlinear terms are irrelevant, i.e. when 



K < Kf < Kg = Kx ■ 



(54) 



As the temperature decreases, the height surface be- 
comes more rigid. There are several possible scenarios 
for a transition from a high-temperature (rough) to a 
low-temperature (smooth) phase. One is the fluctuation- 
driven roughening transition, involving the system freez- 
ing into a smooth phase due to the last term in (39) when 
K approaches Kf from below. 

Based on the evidence accumulated in our numerical 
simulations we conclude that such a behavior does not 
occur. We observe that the measured stiffness k obeys 
the condition (54) as long as the height surface fluctu- 
ates appreciably. Hence, in particidar, the coupling (40) 
to the gate voltage is irrelevant in the correlated phase of 
the model with the long range interaction (2). This con- 
clusion generalizes the observation^* that the coupling to 
external field is irrelevant in the AIAFM (for sufficiently 
small field). 

Lowering the temperature further produces a sharp 
freezing into a ground state in which the height surface 
does not fluctuate. While the freezing is accompanied by 
a steep increase of k to immeasurably high values, the 
latter takes place at T w Tc. Due to narrow temperature 
interval in which this increase occurs, it is difficult to at- 
tribute it to the effect of height fluctuations. Instead, the 
observed behavior is consistent with an ordinary freezing 
by a first or second order phase transition (as discussed 
in Section 111). 

Above the freezing temperature, in the correlated 
phase the density and height fields effectively decouple, 
since the coupling A is irrelevant. Indeed, the inequality 
(48) cannot be satisfied in the correlated phase, Eq. (54). 
This justifies a posteriori our phenomenological approach 
of including both the height and the density in the par- 
tition function (37) . Similar arguments also explain why 
the negative correction (52) to the stiffness does not cause 
an instability. The latter would require k> k\, the con- 
dition forbidden in the correlated phase by Eq. (54). 

Finally, the absence of the dislocation-unbinding phase 
transition between the correlated fluid and the high tem- 
perature disordered phase (studied in Section Vll above) 
can also be understood based on the condition (54). Since 
K < Kf < kkt, the lower bound (36) is never reached. 
[Similar argument was used in Ref. 34 for the AIAFM in 
zero or small external field]. 



IX. DYNAMICS OF THE HEIGHT FIELD 

The weakly coupled fluctuations of height and density 

in the correlated phase lead to interesting dynamical ef- 
fects. Employing the Langevin dynamics, here we dis- 
cuss how the interplay between the h and n fields affects 
the dynamical properties. We find the corrections to the 
conductivity and compressibility due to the height fluc- 
tuations, perturbative in the height-density coupling A. 

The non- conserving Langevin dynamics for the height 
field has the form 



dth{r,t) = -r]T^-^+Cr,t, 



(55) 



with given by Eq. (38), and ^ the stochastic force, 

{^r,t^r',t') = {e)S{t-t')5{T-r'). (56) 

Since the temperature T is included in it multiplies 

the kinetic coefficient in Eq. (55), so that the form of 
the fluctuation-dissipation relationship is preserved: 



(57) 



The conserving Langevin dynamics for the charge density 
n is deflned using the continuity equation 



dtn + Vrj = 0, j = -(TTVr — -hj^, 

on 

with the fluctuating extraneous part obeying 



(58) 



(i^(r,t)j^(r',i')) = ((j'')') 5,, 5{t-t') Sir -r'). (59) 

The conductivity a is related to variance by the 
Nyquist formula 



((jY)=2<7T. 
Eqs. (55, 58) and (39, 41, 42) yield 

27rg 



(60) 



dth{r,t) = j^Ti^KV^h+'^ sin^ 

iriXn + f) . TTh\ . 
+ sin— j +C(r,t), (61) 



dtn{r, t) = aTVl (^j dPr' Ur-r'Ur' + 

-Vr(f+j^), (62) 

where we write the stochastic contribution to the current 
due to the height field fiuctuations in the form 



j''(r, t) = -Ao-TVrCos 



Trhr 



(63) 



The height field fluctuations provide an additional con- 
tribution Sa'* to the conductivity, which is determined by 
the fluctuation-dissipation theorem: 

{j';{-u;,-k)j':{w,k)) = 2Ja;:,(a;,k)T. (64) 
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The contribution to the eonductivity that arises from the 
height field fiuctuations is purely longitudinal: 



a' {-uj, -k) cr' {uj, k)k^ku F{uj,'k) . 



(65) 

Here C7'(u;,k) is the longitudinal part of the total con- 
ductivity a, and F(a;,k) is a Fourier transform of the 
correlator 



Fiv,t) = 



cos ; COS -— 





1 r 7r2 



(/lr,t — hofi)' 



(66) 



Averaging in (66) with respect to the Gaussian free en- 
ergy, A = 5 = in Eq. (61), yields 



(/lr,t — hofi)' 



A\z.r—iujt I 



-/ 

TTK Jo 



1/a , 



k 



• (67) 



Wc arc interested in the asymptotic behavior of the 
height correlator at large separation r ^ a. From the 
small k expansion under the integral (67) we have 



2 /r^ 
1 _ e-'f-fc 1*1 Jo(fcr) ~ f — + r]K\t 



(68) 



Therefore the height correlator (67) asymptotic behavior 
at r a is consistent with Eq. (35): 



(69) 



At large r, t the correlation function (66), as well as the 
current j'' correlation (64), (65), are thus of a power law 
form: 



F(r,i) = 



2 VrV4 + ?7K|i| 



Ik* Ik 



(70) 



With K < K* in the correlated phase, the function F(r, t) 
rapidly decays at large spatial and temporal separations. 
We conclude from Eqs. (64), (65) and (70) that in this 
case the current j'*(r, t) space and time correlations have 
short memory and are local. Hence one may treat j''(r, t) 
as an additional source of the Johnson-Nyquist noise, and 
our approach based on the fluctuation-dissipation theo- 
rem (64) is a posteriori justified. 

Finally we consider the correction to the compressibil- 
ity of the charged system due to the height fluctuations. 
It can be obtained directly from the statistical averaging 
with respect to the canonical distribution. The compress- 
ibility u, defined as 



(<l>_k^>k) , with $ = 



Sn 



(71) 



acquires the perturbative correction as a result of the 
height-density coupling (42). The uniform {k = 0) part 
Su is obtained from 

^ ^ Uk^o + SUk=o , (72) 
with the height field-induced density- density interaction 



cos — ;— cos — ; — ) . (73) 
b b /o 



The correlator (73) is given by Eq. (70) with t = 0. Thus 
the compressibility correction due to fluctuations 



6iy = -X^ [C/k=o]"^ j (fvF{r,t = 0) . 



(74) 



This correction is finite in the correlated phase. Indeed, 
the correlator in (73) is local when k < 2k* which is 
consistent with the condition (54), and, thus, the integral 
in Eq. (74) is infrared-convergent. 



X. SUMMARY 

In the present work we analyze the problem of charge 
ordering and dynamics of classical electrons on a 2d tri- 
angular array, by relating it to the better studied prob- 
lem of the triangular Ising antiferromagnet. The phase 
diagram (Fig. 3) is found to be more rich than in the 
latter problem due to long-range electron interactions in- 
terplay with the geometrical frustration. At low temper- 
atures, the electron system freezes into commensurate or 
disordered ground states in a continuous range of densi- 
ties. We demonstrate that transport can be employed to 
study charge ordering, with the singularities of conduc- 
tivity marking the phase transitions. 

At intermediate temperatures, we idenitify the topo- 
logical fluid phase characterized by strong electron cor- 
relations. The fluid state is described in terms of a non- 
local order parameter, the height field. We find that 
the short range correlations in the charge system correp- 
sond to Gaussian fluctuations of the height fleld in the 
presence of topological defects. We determine numeri- 
cally the effective stiffness of the height field, and rule 
out the Berezinskii-Kosterlits-Thouless phase transition, 
concluding that the topological defects remain unbound 
above freezing. In this phase the electron transport, re- 
sulting from short-range correlations enforcing local con- 
tinuity of the height field, takes the form of "free" charge 
dynamics corresponding to height field fluctuations. We 
found that both the topological defects and the height 
fluctuations contribute to transport, their relative con- 
tribution depending on the concentration of the defects. 

The ground state properties have been studied in this 
work for the simplest fractions n = 1/3, 2/3, and n = 
1/2. In transport, freezing is manifest by singularities 
in zero bias conductivity which drops to zero in the or- 
dered phases. The sensitivity of transport to ordering of 
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the electron system makes it a useful probe of 2d charge 
states. 

While the stuation at n — 1/3, 2/3, 1/2 is fairly 
conventional, the precise nature of ordering for generic 
density remains unclear. Here several different scenar- 
ios can be anticipated. One is a devil's staircase of in- 
compressible states at any rational n. In the relatively 
small system studied using MC dynamics, we have in- 
deed found freezing into higher order fractions with rela- 
tively small denominator, such as the "striped" n = 3/7 
state shown in Fig. 6 (D). Similarly, other states with 
rational n = p/q could become incompressible at low 
temerature, with decreasing with the increase of the 
denominator q, forming the devil's staircase. Another 
possibility is the appearance, besides the simple rational 
fractions, of a family of disordered ground states. One 
documented example of such a behavior is the prerough- 
ening phenomenon'^^ which describes a continuous phase 
transition into a disordered flat phase of a 3d crystal 
surface, characterized by the height surface with a dis- 
ordered array of steps. The tunability of the form of 
interaction in the quantum dot arrays should allow to 
explore these, and other complex states. 
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APPENDIX A: ASSIGNING HEIGHT FIELD TO 
THE ARRAY 

Here we describe the procedure by which the height val- 
ues are assigned to the sites of the array. It should be 
noted that the height is well defined only at T = 0, in 
the absence of topological defects, while at finite tem- 
perature it is defined only locally, in the regions away 
from the defects. As we found in Section V, the defects 
are always present in the correlated phase, since there is 
no defect binding phase transition. This complication, 
however, is inessential, since in the number of defects is 
exponentially small at low enough temperature, and is 
often zero for the finite patch used in our simulation. We 
found that defining the height locally in a system with a 
small number of defects does not lead to any significant 
errors in the statistics of fiuctuations and, in particular, 
does not affect the numerical value of the effective stiff- 
ness AC. 




FIG. 10: (Color online) Assigning auxiliary 3d coordinates to 
the sites of the array. 

We work with an N x N rhombic patch of the 2d trian- 
gular array, using the coordinate system aligned with the 
array, as shown in Fig. 10. The sites are labeled by integer 
coordinates i,j — 0, — 1. We place the origin 

in the upper left corner of the patch. The first compo- 
nent, i, is the site number counted along the horizontal 
axis (the patch upper edge). The numbers i increase as 
we go from left to right. The second component, j, is the 
site number along the axis which points downward and to 
the right at the angle 7r/3 with the horizontal axis. The 
numbers j increase as we go from the origin down and to 
the right, with j = N ~ 1 a,t the lower edge of the patch. 
The conventional Cartesian coordinates — (x^mym) 
of the site {im,jm) are given by 

— • ■ ^ — • ^ ( ^^^ 

•^m — ^ ~^ Jm i Um — Jm ^ ^ • \A1) 

To calculate the height h{rm) for each point r„j = 
{xrmym) we first assign the auxiliary 3d coordinates 

R(r„) = ^i?,(r„Oe,. (A2) 

s 

Here the unit vectors e^, s = 1, 2, 3, are the basis vectors 
of the auxiliary 3d simple cubic lattice. 

To obtain the 3d coordinates R(i,j) for a particular 
charge configuration, we start from the upper left corner, 
R(i = 0,j = 0) = 0, and assign the coordinates by the 
sequence of steps. We shall call the sites (i,j) and («', j') 
connected, if the charges at these sites are opposite, i.e., 
the bond between these sites is not frustrated (and hence 
is drawn in Fig. 10). 

To assign the height, we move from point to point in 
each row, i = 0, A^ — 1, with j fixed, repeating it for all 
j = 0, A — 1. Suppose that the heights of all the sites 
with i' < i and j' < j are already defined. The heights of 
the remaining one, two, or three sites that are connected 
to (z, j) are defined by the following rules (Fig. 10): 

(i) If the site + is connected to then Ri{i + 
= Ri{i,j) + b; 

(ii) If the site j + 1) is connected to then 
R2{i,j + l) = R2{i,j)-b; 

(ni) If the site {i — l,j + 1) is connected to (i, j), then 
i?3(i -IJ + 1) = RsiiJ) + b, with b given by Eq. (24). 
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Once the 3d coordinates R(rm) are defined, the height 
at each site r„ is obtained by projecting onto the (111) 
crystal direction: 

3 

h{Vm) = J2 ^«(''™) • (A3) 

.s=l 

We note that for this procedure to work, each site must 
be connected, in the above sense, to at least one neigh- 
bor. Thus the height assignment may be impossible for 
charge configurations with very fow density of unfrus- 



trated bonds. However, at low enough temperatures, in 
the correlated fluid phase, at n ~ 1/2, we do not en- 
counter such configurations in our simulation. 

Also, in the presence of defects, the above recipe as- 
signs heights unambiguosly, albeit in a somewhat ad hoc 
way. We do not explicitly exclude the configurations with 
defects from the height statistics analysis. Instead, we es- 
timate their relative contribution (Sec. VI) and find it to 
be small enough for the results (the stiff'ness value) to be 
affected. 
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